Field of the Invention
The present invention relates to a solution method and a solution apparatus for solving a system of linear equations, and more particularly to a technique of solving an underdetermined system of linear equations in which the number of equations is smaller than the number of variables.
Description of the Related Art
In a medical instrument used for X ray computed tomography (CT), magnetic resonance imaging (MRI), or the like, a three-dimensional structure of the inside of a human body is reconstructed from an observation image projected onto a two-dimensional plane. When a three-dimensional structure to be estimated is set as an unknown vector x, an observation process such as X ray projection is set as a matrix A, and an observation signal is set as b, the observation process can be formulated as follows:Ax=b  (1)
Here, A may also be referred to as a coefficient matrix, and b may also be referred to as a constant vector or a right side vector.
When solving Equation (1), the number of observations is typically increased in order to improve the estimation precision of the three-dimensional structure. Hereafter, a dimension of a signal x to be estimated is set as N, a dimension of the observation signal b is set as M, and it is assumed that M<N. In Equation (1), the number M of equations is smaller than the dimension N of the variable, and therefore Equation (1) is referred to as an underdetermined system of linear equations.
An underdetermined system of linear equations has an infinite number of solutions, and therefore a solution is typically determined uniquely by applying a technique known as regularization.
In a recently proposed technique, an underdetermined system of linear equations is solved with a high degree of precision by assuming that the signal to be estimated is sparse (in other words, assuming that a majority of elements of the vector x are zero). This technique is known as sparse regularization, and is currently attracting attention as an effective technique for removing noise from images, reconstructing medical images, and so on, for example.
Algorithms of sparse regularization can be broadly divided into two types, namely a “greedy algorithm” and “L1 regularization”. A method exhibiting the most favorable performance with a greedy algorithm is subspace pursuit (abbreviated at SP hereafter), disclosed in Non-patent Literature (NPL) 1, and a method exhibiting the most favorable performance with L1 regularization is FISTA, disclosed in NPL 2.
A flow of processing performed in SP, with which a more precise solution can be found at a much higher speed than with FISTA, will now be described using FIG. 3.
In step S301, values of the coefficient matrix A and the constant vector b are set. Further, the number of non-zero elements K is set. The value of K is determined by a user, and K may be set at any value that satisfies K≦M/2.
Below, a set of element numbers of the variable x is denoted as Ω, and a universal set is denoted as Ωall:={1, 2, . . . , N}.
In step S302, the variable x(0) is initialized. Further, a set Ω(0) of element numbers of non-zero element candidates is initialized as an empty set. The number of iterations t is set at 1. Note that the characters in superscript parentheses express the number of iterations.
In step S303, a residual r(t−1):=b−Ax(t−1) is calculated using the variable x(t−1).
In step S304, a pseudo-error ε(t−1):=ATr(t−1) is calculated, and K elements are extracted in order from the element having the pseudo-error with the largest absolute value. A set of element numbers of the K extracted elements is set as ΔΩ(t).
In step S305, a set Ωp(t)=Ω(t−1)∪ΔΩ(t) is determined, and a projection operator P(t):=Ωall→Ωp(t) is generated. The number of elements in Ωp(t) is no smaller than K and no larger than M. For example, when Ωall:={1, 2, . . . , 5} and Ωp(t):={1, 3, 5}, the projection operator P(t) takes the form of a following matrix:
                              P                      (            t            )                          =                  (                                                    1                                            0                                            0                                            0                                            0                                                                    0                                            0                                            1                                            0                                            0                                                                    0                                            0                                            0                                            0                                            1                                              )                                    (        2        )            
In step S306, a subproblem expressed by a following equation is constructed and solved.(A(P(t))T)z(t)=b  (3)
Here, z(t):=p(t) x(t) is a variable of the subproblem, a dimension of which is no smaller than K and no larger than M.
When the dimension is M, Equation (3) is regular such that a solution is determined uniquely, and therefore a normal solution method for a system of equations, for example a direct method such as the Jacobi method, the Gauss-Seidel method, or the SOR method, or an iterative solution method such as the conjugate gradient method, is used. When the dimension is smaller than M, on the other hand, this indicates an overdetermined system, and therefore the subproblem is considered as a problem of the method of least squares. Accordingly, a following equation, for example, is solved:(P(t)ATA(P(t))T)z(t)=P(t)ATb  (4)
In step S307, a value of x(t) is calculated from a following equation using a solution z(t) obtained in step S306.x(t)=(P(t))Tz(t)  (5)
In step S308, K elements are extracted from x(t) in order from the element having the largest absolute value, and a set of element numbers of the K extracted elements is set as Ω(t). Ω(t) is used in step S305 of the next iterative calculation.
Steps S303 to S308 are then repeated until a predetermined convergence condition is satisfied.
[Non-patent Literature 1] Wei Dai and Olgica Milenkovic: “Subspace Pursuit for Compressive Sensing Signal Reconstruction,” IEEE Transactions on Information Theory, VOL. 55, NO. 5, pp. 2230-2249, MAY 2009.
[Non-patent Literature 2] Amir Beck and Marc Teboulle: “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems,” SIAM Journal on Imaging Sciences, Vol. 2, No. 1, pp. 183-202, 2009.